Turbulent thermalization of the Quark Gluon Plasma 
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Classical-statistical lattice gauge theory simulations are employed to demonstrate the existence 
of a nonthermal fixed point in the space-time evolution of heavy ion collisions at ultrarelativistic 
energies. After an initial transient regime dominated by plasma instabilities and free streaming, the 
ensuing overpopulated non-Abelian plasma exhibits the universal self-similar dynamics characteristic 
of wave turbulence observed in a large variety of physical systems across different energy scales. 
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The past decades have revealed remarkable properties 
of quark-gluon matter formed in heavy-ion collision ex- 
periments at RHIC and the LHC, in particular its ap- 
parent behavior as a nearly perfect fluid. Understanding 
how these features emerge ab initio in a heavy-ion col- 
lision, and to what extent a thermalized Quark Gluon 
Plasma is achieved, is an outstanding problem in theo- 
retical physics pp. 

Some of the out-of-equilibrium properties of non- 
Abelian plasmas may be shared by a large class of differ- 
ent quantum many-body systems on the most diverse en- 
ergy scales. Universal behavior far from equilibrium has 
been predicted for systems ranging from early-universe 
inflaton dynamics to table-top experiments with cold 
atoms. The universality is based on the existence of non- 
thermal fixed points, which are attractor solutions with 
self-similar scaling behavior. In these examples, they are 
associated with the phenomenon of wave turbulence [2j- 
5 . While such a behavior in weakly coupled non-Abelian 
plasmas without expansion is rather well understood [5]- 
18], the anisotropically expanding plasma encountered in 
high-energy heavy-ion collisions is much less so. 

In recent years there has been significant progress in a 
first-principles understanding of expanding non-Abelian 
plasmas in two limiting cases. One of these is the 
study of the strong-coupling limit using gauge-string 
dualities in supersymmetric Yang-Mills theories. The 
results in this case indicate the important role of 
anisotropies in the longitudinally expanding system even 
at the transition to the hydrodynamic regime 9 . The 
other case that is amenable to ab initio calculations is 
Quantum Chromodynamics in the weak-coupling limit 
a s < 1. The colliding nuclei in this limit are described 
in the Color Glass Condensate framework [10J . The 
dynamics of the non-equilibrium Glasma created in 
such a collision is that of highly occupied gluon fields 
with typical momentum Q. Since the characteristic 
occupancies ~ l/a s (Q) are large, the gauge fields are 
strongly correlated even for small gauge coupling and 
the Glasma exhibits classical-statistical dynamics. 



In this work, we address the problem of thermaliza- 
tion of the Quark Gluon Plasma in the weak-coupling 
limit. Within this framework, a number of scenarios of 
how thermalization proceeds have been developed [T2T 
I15j . However, the available parametric estimates allow 
for different types of solutions of the employed kinetic 
equations. Since classical-statistical lattice gauge theory 
and kinetic theory have an overlapping range of validity 
[15] , together they can be used to obtain an unambigu- 
ous description of the dynamics starting from the non- 
perturbative early stages with occupancies ~ 1 /a s to the 
late-time quantum regime. 

With this goal in mind, we perform classical-statistical 
lattice simulations of SU(2) gauge theory in 3 + 1 dimen- 
sions with longitudinal expansion. We discover that after 
a transient regime the dynamics becomes independent 
of the details of the initial conditions and leads to the 
same attractor solution. In this regime, the plasma ex- 
hibits self-similar behavior characteristic of wave turbu- 
lence and we compute the corresponding universal scal- 
ing exponents and scaling functions. A striking conclu- 
sion is that while the physics of plasma instabilities and 
free streaming is found to play an important role for the 
dynamics at early times, they do not govern the uni- 
versal turbulent regime. Instead, the attractor solution 
we find shows the same characteristic time-dependence 
as the "bottom-up" thermalization scenario [T5] to high 
accuracy in the range of applicability of our simulations. 

The expansion of the plasma in the longitudinal di- 
rection x 3 at time x° is conveniently discussed in terms 
of the proper time t = W (x ) 2 — (x 3 ) 2 and rapidity 
i] = atanh(x 3 /x°), and we will denote the transverse 
coordinates by xt = (x^x 2 ). For our simulations we 
employ the Hamiltonian lattice formulation for SU(2) 
gauge theory in Fock-Schwinger gauge (A t — 0) fol- 
lowing Refs. [TTll20| . In this framework, gauge in- 
variant observables can be calculated from the stan- 
dard lattice plaquettes U^ v {x), which depend on space- 
time variables x^ for fi,v = (t,x\,X2,T)). The plaque- 
ttes are related to the non-Abelian field-strength ten- 
sor F fJ _„(x) = d ll A u - d v A^ + ig[A^,A u ] in terms of 
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FIG. 1. Ratio of longitudinal to transverse pressure as a 
function of time for different initial anisotropies £o and fixed 
initial occupancy no = 1. The inset shows the same quantity 
for different initial occupancies no and fixed initial anisotropy 
£o = 1 along with the free streaming (dashed) curve. 
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FIG. 2. Time evolution of the characteristic longitu- 
dinal momentum scale for different initial anisotropies £o- 
The inset shows the scaling exponent extracted for differ- 
ent initial conditions and lattice discretizations. The average 
27 = 0.67 ± 0.07 is indicated by gray lines. 



U(j,v(x) — exp^ga^a^F^u) for sufficiently small lattice 
spacings a M , where g 2 = 47ra s . 

In order to extract gauge dependent quantities, we 
use the residual gauge freedom to perform time inde- 
pendent gauge transformations such that the generalized 
Coulomb gauge condition d^Ai + t~ 2 d n A r) = (i=l,2) is 
imposed at each read-out time. In particular, we use 
this procedure to extract the gluon distribution func- 
tion f(pT,p z ,t), which describes the gluon occupancy 
per momentum mode, averaged over spin and color. Here 
we identify the longitudinal momentum at mid-rapidity 
as p z = v/t, where v is the Fourier conjugate with re- 
spect to relative rapidity, and pt denotes transverse mo- 
mentum. Our initial Gaussian probability distribution 
of gauge field configurations on the lattice can also be 
characterized in terms of the initial gluon distribution 
function 

f(PT,Pz,t ) = ^ Q (Q-\Ipt + , (1) 

which describes the overpopulation of gluon modes up 
to the momentum Q. The parameter hq is used to vary 
the initial occupancy, while £0 parametrizes the initial 
anisotropy in momentum space. The initial time is cho- 
sen as Qto = 100 to minimize discretization errors, while 
accessing sufficiently late times t ^> t in order to observe 
universal aspects of the evolution. We explicitly verified 
that varying the initial time in a range up to Qto = 1000 
only affects the transient regime but does not change the 
universal properties at later times. In general, earlier ini- 
tial times strengthen the role of free streaming, while for 
larger Qto the transient regime starts to resemble static- 
box simulation results [SHE]- 



A central gauge-invariant quantity is the stress-energy 
tensor, which is expressed in terms of field strengths as 
T^ v = -g va F* 5 F aS + g^F"i s F l5 /A with the non-trivial 
metric g^ u — diag(l, — 1, — 1, — t 2 ). The transverse and 
longitudinal pressures are 

PT{t) = -\(T x x {t)+Ty{t)) ,P L (t) = -(l*(t)), (2) 

where the brackets denote Monte Carlo averaging with 
respect to the initial probability distribution. In Fig.[l]we 
show the ratio Pt / Pt as a function of time. The curves in 
the main graph are for different initial anisotropies £0 an d 
fixed initial occupancy no = 1. Starting from an isotropic 
initial distribution (£0 = l)j the system is seen to become 
more and more anisotropic with time as a consequence 
of the longitudinal expansion. Indeed, the early-time be- 
havior is governed by free streaming, whereas at later 
times the anisotropy of the system increases more slowly 
as a consequence of interactions. This is elaborated fur- 
ther in the inset where the free streaming (dashed) curve 
is shown for comparison, along with results for different 
initial occupancies no and fixed £0 = 1- For strong ini- 
tial anisotropy, such as for £0 = 4 and 6, there is a short 
transient regime where Pl/Pt increases due to plasma 
instabilities. Most remarkably, this does not affect the 
evolution at later times: After the transient regime all 
curves show very similar scaling with time, irrespective 
of the choice of initial conditions. 

In order to analyze this behavior in more detail, we in- 
troduce the transverse and longitudinal hard momentum 
scales At and A^. These are gauge- invariant quantities 
and are computed from longitudinal and transverse pro- 
jections of the square of the covariant derivative of field 
strengths divided by the energy density [Ej- 
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FIG. 3. Rescaled moments of the distribution function as a 
function of the rescaled longitudinal momentum. The rescaled 
data for different times all collapse onto a single curve. 

A universal behavior of these hard scales can be char- 
acterized in terms of scaling exponents 7 and /3 as 

A 2 L (t)~(Qt)- 2 \ A 2 T (t)~(Qt)- 2 ?. (3) 

In Fig. |2j the time evolution of the rescaled K 2 L for differ- 
ent initial £0 and fixed no = 1 is shown. After the tran- 
sient regime, for Qt > 650 one clearly observes a univer- 
sal power-law behavior. Comparison to the dashed curve 
~ (Qt)~ 2 / 3 indicates an approximate value of 7 ~ 1/3. 
The inset of Fig. [2] gives the fitted scaling exponent of 
K 2 L as a function of time for a set of four different initial 
conditions in the range £0 = 1 - 6 and uq = 0.25 - 2. 
To check for a possible further dependence on lattice dis- 
cretizations, we also display results from the evolution 
f° r Co = no = 1 using four different N T x N v lattices 
in the range N T = 256 - 512, N v = 1024 - 4096 with 
Qa T = 0.5-1 and a v = (0.625 - 2.5) • 10~ 3 . By aver- 
aging over the data points for Qt > 650 we obtain the 
estimate 27 = 0.67 ± 0.07. 

In contrast, we find that in this regime At stays ap- 
proximately constant in time, which corresponds to j3 ~ 
0. Indeed, the exponent extracted from a fit to our data 
for A^ approaches zero monotonically with \f3\ < 0.04 
for Qt > 650. In general, we observe that the errors are 
not dominated by the finite lattice but the remaining de- 
pendence on the initial conditions for the available finite 
times. We have also verified that the Debye scale, which 
is relevant also for the physics of plasma instabilities, is 
resolved for all simulation times on the large lattices. 

In terms of the gluon distribution function, a self- 
similar evolution has to fulfill 

f( P T,p z ,t)=t a f s (tPp T ,t~<p z ), (4) 

where fs denotes a stationary distribution independent 
of time. The scaling exponents a, /3 and 7 are universal, 
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FIG. 4. Dependence of the gluon distribution function on 
transverse momentum pr for p z = at different times. The 
inset shows the self-similar evolution for the rescaled second 
moment of the distribution as a function of pr for p z — 0. 

as is the stationary scaling function up to amplitude nor- 
malizations. The scaling exponents j3 and 7 agree with 
the previous definition ([3]), as can be seen by evaluating 
the perturbative expression for the hard scales - consid- 
ering only the Abelian parts of the field strength tensor 
- as 

A 2 / t j _ I d2 PT I d Pz Pt. z uj p f{pr,Pz,t) 
T ' L J d 2 p T J dp z uj p f(pT,Pz,t) 

Here ui p ~ pt is the relativistic quasi-particle energy 
in the scaling regime, where the characteristic trans- 
verse momenta are much larger than the longitudinal 
momenta. 

In order to investigate the emergence of a self-similar 
scaling solution in our simulations, we study the time 
evolution of the gluon distribution function. Fig.[3]shows 
the rescaled zeroth and second moment of the gluon dis- 
tribution as a function of fjp z for pt — Q, where we 
define t r = Qt/10 3 . The rescaled data for different times 
in the range from Qt = 750 to 4000 is seen to collapse 
onto a single curve. This is a striking manifestation of 
self-similarity. Fig. [4] displays the distribution, now as a 
function of transverse momentum for p z — at different 
times. The spectrum can be characterized as a thermal- 
like 1 Ipt power-law over a large range of transverse mo- 
menta pt < At with a rapid fall-off for pt > At- In- 
deed, the position of At is seen to remain approximately 
constant in time and the spectrum exhibits a self-similar 
evolution with a decreasing amplitude 

n Ha rd(*) = fipT = Q,Pz = 0)*)- (6) 

In the inset of Fig. |1J the rescaled second moment of 
the distribution, ~ t~ a p T f(pT,p z = 0,t), is shown on a 
linear scale. The results at different times again nicely 
collapse onto a single curve to good accuracy. 
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In order to determine the scaling exponent a, we have 
employed a simultaneous least squares fit to self-similar 
behavior of the zeroth to fourth moment of the distribu- 
tion function. Here it is advantageous to estimate the 
combination a — 7. Our data shows that a — 7 = — 1 — A 
with a monotonically decreasing A < 0.13 for Qt > 650. 
Taking into account our above estimate for 7, this leads 
to a = —0.67 ± 0.04 — A. Comparing the statistical un- 
certainty to the systematic error A, one observes that the 
dominant source of error is again the fact that the scal- 
ing limit is only approximately realized for the available 
finite times. 

We will now provide a simple explanation of the char- 
acteristic scaling behavior observed in our simulations. 
In the framework of kinetic theory, the time evolution 
of the gluon distribution is described by a Boltzmann 
equation of the form 



Pz 



f(pr,Pz,t) = C[p T ,Pz,t;f] , (7) 



with a generic collision term C[pT,p Zl t; f] for n -H- m 
scattering processes. We follow here the standard tur- 
bulence analysis [2] for self-similar evolutions of the type 
given by Q. In this analysis, the scaling behavior of the 
collision integral 



C\p T ,p z ,t;f}=^C[^ PT ,r Pz ;f s ] 



(8) 



is described in terms of the exponent /1 = /i(a,/3,7), 
whose precise form depends on the underlying interac- 
tion. Substituting this form into the Boltzmann equation 
|7]), one finds that the fixed point solution fs{PT,Pz) sat- 
isfies the relation 

afs(pr,Pz) + f3p T d PT fs(pT,Pz) 

+ (7- l)Pzd p Js{pT,Pz) = C[p T ,p z ;fs] , (9) 



with the scaling condition, 

a - 1 = /x(a,/3,7) 



(10) 



As previously argued in the "bottom-up" scenario 
(henceforth labeled BMSS) [12], the interaction of hard 
excitations is dominated by elastic scattering with small 
momentum transfer, as long as their occupancies are 
large (njiard ^ !)• The dominant effect of these collisions 
on the particle distribution f(pT,p Zl t) is to broaden the 
longitudinal momentum distribution by multiple incoher- 
ent small-angle scatterings. This may be characterized by 
a collision integral of the Fokker-Planck type, 

C (olast) [PT,p z ;/] = q d 2 p J(p T ,Pz,t) . (11) 

The momentum diffusion parameter q in this expression 
is parametrically given by 



2 / d2 PT f dPz , 2 , ,\ 

■ 1 ' f (PT,Pz,t) 



2tt 



(12) 
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FIG. 5. Evolution in the occupancy-anisotropy plane. In- 
dicated are the attractor solutions proposed in (BMSS) [12] . 
(KM) 14 and (BGLMV) [T5], along with the simulations re- 
sults for different initial conditions (in blue). 



in the limit of high occupancies. This approximation is 
supposed to describe the dominant physics relevant for 
the overall scaling with time, which enters the scaling 
relation (10 1 considered in the following. However, it 



does not have to be an accurate approximation for the 
solution of the fixed-point equation Q , which in general 
requires more detailed information about the momentum 
dependence of the collision integral. 



For the collision integral in ( 11 ), the self-similar distri- 
bution Q leads to fi = 3a - 2/3 + 7. Together with ( [Io| , 
one obtains the scaling relation 2a — 2/3 + 7 + 1 = 0. Since 



elastic scattering (11) is particle number conserving, a 



further scaling relation is obtained from integrating the 
distribution function over px and rapidity wave numbers 
v = tp z . By use of the scaling form Q, this constraint 
leads to the scaling relation a — 2(3 — 7 + 1 = 0. A cor- 
responding scaling relation can be extracted from energy 
conservation. Taking into account that the mode energy 
lj p ~ pt in the anisotropic scaling limit, this yields the 
condition a — 3/3 — 7 + 1 = 0. The three scaling relations 
from small-angle elastic scattering, particle number and 
energy conservation, as incorporated in the BMSS kinetic 
approach give 



-2/3 



= 0, 7 = 1/3, 



(13) 



which is in excellent agreement with our lattice simula- 
tion results. 

Besides the BMSS scenario, alternative thermalization 
scenarios at weak coupling with different attractor so- 
lutions have been proposed in the literature. In the 
(KM) scenario [2], plasma instabilities play a key role 
for the entire evolution, while in the (BGLMV) scenario 
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[T5] elastic scattering and transient Bose-Einstein con- 
densation is argued to generate an attractor with fixed 
anisotropy parameter S s . Our findings are compactly 
summarized in Fig. [5j describing the space-time evolu- 
tion in the occupancy-anisotropy plane. The horizontal 
axis shows the occupancy nnard and the vertical axis the 
momentum-space anisotropy in terms of the typical lon- 
gitudinal and transverse momenta At,l- The gray lines 
indicate the different attractor scenarios while the blue 
lines show our simulation results for different initial con- 
ditions. One clearly observes the attractor property, in 
very good agreement with the analytical discussion of the 
BMSS kinetic equation in the high-occupancy regime. 

In the quantum regime, where occupation numbers be- 
come of order unity, the classical-statistical framework 
can no longer be applied. In the BMSS scenario, in- 
elastic processes begin to play a significant role in this 
regime and eventually lead to thermalization. Here it is 
also conceivable that effects such as plasma instabilities 
play a role for the thermalization process at this stage of 
the evolution [151 HI]. 

We note that for realistic values of a s at RHIC and 
LHC energies, the dynamics takes place on much shorter 
time scales where details of the initial conditions may 
be essential [2T| [22] . Since the time for the approach to 
the universal attractor is controlled by the initial condi- 
tions, it will be interesting to explore how and on what 
time scales the results of [2T] converge to the universal 
dynamics discussed in this work. 

Finally, we observe that in anisotropic hydrodynamic 
descriptions [35J of heavy ion collision data, the essen- 
tial ingredient is not isotropization or thermalization but 
rather the existence of an anisotropic attractor solution. 
If the concepts developed here are applicable in heavy ion 
collisions at RHIC and LHC energies, non-equilibrium 
wave turbulence and attendant phenomena, may be play- 
ing a larger role in the dynamics of the Quark Gluon 
plasma than earlier conceived. 
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